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To analyze possible generalizations of reaction-diffusion schemes for the case of subdiffusion we 
discuss a simple monomolecular conversion A — » B. We derive the corresponding kinetic equations 
for local A and B concentrations. Their form is rather unusual: The parameters of reaction influence 
the diffusion term in the equation for a component A, a consequence of the nonmarkovian nature of 
subdiffusion. The equation for a product contains a term which depends on the concentration of A 
at all previous times. Our discussion shows that reaction-subdiffusion equations may not resemble 
the corresponding reaction-diffusion ones and are not obtained by a trivial change of the diffusion 
operator for a subdiffusion one. 
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Many phenomena in situations out of equilibrium can 
be described using a picture based on reaction processes. 
Apart from chemical reactions, the examples are exciton 
quenching, recombination of charge carriers or radiation 
defects in solids, predator-pray relationships in ecology 
etc. Reactions in homogeneous media are often described 
by formal kinetic schemes. Thus, the concentrations 
Ci(t) of the components follow the first-order differen- 
tial equations dCi{t)/dt = fi{C\{t), CV(i)} where the 
reaction terms typically have a form fi{C\, Cjy} = 
±KiC^ 1 C% 2 ...C™ N with the powers rij depending on the 
stoichiometry of the reaction and Ki denoting the cor- 
responding reaction rates. In inhomogeneous situations 
(layered systems, fronts, etc.) the mesoscopic approach 
based on reaction-diffusion equations for the position- 
dependent concentrations Ci (r, t) is often the appropriate 
way of description. In case of normal diffusion such equa- 
tions are obtained by adding a diffusive term to classical 
reaction schemes and have the form 



dCj(r,t) 

at 



KiAdiv, t) + fi 



(1) 



with fi = /j{Ci(r, f), Cjy(r,t)} and Ki being the dif- 
fusivity of the component i. This approach is applicable 
whenever characteristic scales of spatial inhomogeneities 
are much larger than the typical interparticle distances 
and particles' mean free paths (see e.g.0). As we pro- 
ceed to show, the possibility to put down such schemes 
is due to the Markovian nature of normal diffusion. 

Nowadays more and more attention is paid to situa- 
tions when the diffusion is anomalous, which are found 
to be abundant 0. One of the most important situa- 
tions here is the case of subdiffusion described within the 
continuous- time random walk (CTRW) scheme 0. In 
this case the subdiffusive nature of motion stems from 
the fact that particles get trapped and have to wait for 
a time t (distributed according to power-law probabil- 
ity density function tp(t) oc t~ 1 ~ a ) until the next step 
can be performed. It was shown that the properties of 



the reaction under such subdiffusion might be vastly dif- 
ferent from those in diffusive systems H 110,0]. The 
microscopic approach of these works aims on the under- 
standing of the situation when the particles performing 
CTRW react on encounter (and don't react as long as 
they do not move). Such situation is pertinent to exci- 
ton quenching in solids, or to transport in ion channels 

In many cases, however, a mesoscopic approach is 
desirable. Such an approach was adopted in the case 
of reactions under superdiffusion due to Levy flights, 
Refs.j^, , where the transport process involved is 
Markovian. The situation with subdiffusion is much more 
subtle due to strongly nonmarkovian character of subdif- 
fusive transport Here two different situations can 
be encountered: either the reaction at small scales is 
also subdiffusion-controlled (like in the models discussed 
above, where particles can only react if a new step is 
made) or it locally follows normal, classical kinetics. This 
last case that we address here is physically relevant since 
it describes reactions in porous media. The situation is of 
extreme importance in hydrology, where the transport in 
catchments is hindered by trapping in stagnant regions 
of the flow, caves and pores on all scales. The transport 
at long times and large scales is adequately described by 
CTRW |l2| . However on small scales reactions take place 
in normal aqueous solutions, so that particles trapped in 
stagnant regions still can react with each other. A meso- 
scopic approach to such a case was adopted in |l-'"{l| within 
a probabilistic scheme, while [lj] tackle this problem by 
using equations of the same form as our Eq. Q where the 
diffusion operator is changed for a subdiffusion one, con- 
taining an additional fractional derivative in time [3, 1 1 5| : 



dd(r, t) 
dt 



(2) 



In this equation a, is the exponent of the anomalous 
diffusion for the component i, oDl~ az is the operator 
of fractional (Riemann-Liouville) derivative, and -?Q, a4 is 



2 



the corresponding anomalous diffusion coefficient. Such 
equations with decoupled transport and reaction term 
were postulated based on the analogy with Eq.0J and 
look quite plausible. In some cases also a reaction term 
has to be modified by applying a fractional derivative as 
suggested by a microscopic model in 

In what follows we derive the reaction-subdiffusion 
equations for the simplest reaction scheme (monomolecu- 
lar conversion A — ► B) corresponding e.g. to radioactive 
decay of isotope A which is introduced into the ground 
water at some place at time t = and is transported 
according to anomalous diffusion. We show that the cor- 
responding equations do not follow a pattern of Eq.(J2J, 
so that the reaction and diffusion terms do not decouple. 

Let us assume for the time being that all properties 
of A and B particles are the same, so that the reaction 
corresponds to a relabeling of A into B taking place at a 
rate n. In what follows we will use one-dimensional no- 
tation, the generalization to higher dimensions is trivial. 
The Eqs.Q for this case read: 



OA 
~dt 



KAA - kA, 



dB 
~dt 



KAB + kA. 



(3) 



with K being the normal diffusion constant. Let 
C(x, t) = A(x, t) + B(x, t) be the sum of concentrations. 
It evolves according to a diffusion equation: 



Both concentrations A and B follow 



A(x,t) 



f C(x,t), B(x,t) = [1 



(4) 



C(x,t). (5) 



To see this apply Laplace-transform to the equations 10 • 
Note that the solution for A(x, t) in the Laplace domain 
is A(x,u) = C(x,u + k). Eq.© reflects the fact, that the 
conversion is independent from the motion of particles, 
so that concentrations of As and of Bs are proportional 
to the overall concentration multiplied by the probability 
for a particle to survive as A or to become B. The same 
argument leads to the conclusion that Eq.© also holds 
in anomalous diffusion, whatever the evolution equation 
for C is. For subdiffusion 



dC(x,t) 
dt 



K a0 D]- a AC{x,t) 



(6) 



so that in Fourier-Laplace domain for the initial condition 
C(x,0) = 5(x) one has C(k,u) = {u + u 1 - a k 2 K a y 1 so 
that, for instance 



A(k,u) 



1 



{u + n) + (u + Ky- a k 2 K a 



(7) 



However, neither the solution of Eq. (J2J) nor the solution 
of the fractional equation with the modified reaction term 
reproduce this result: the simple reaction-subdiffusion 
schemes do not describe the conversion reaction correctly. 



Let us now turn to deriving a correct reaction-diffusion 
equations for our case. Our derivation will follow the 
way of derivation of the generalized master equation for 
CTRW used in ref. based on the ideas of 0- We 
start from a discrete scheme and consider particles occu- 
pying sites of a one-dimensional lattice. The generalized 
reaction (sub-)diffusion equations follows from the bal- 
ance conditions for particle numbers. A balance condi- 
tion for the mean number Ai of particles A on site i of 
the system reads 



dMt) 
dt 



= i+(t)-i-(t)-KMt) 



(8) 



where JT (t) is the loss per unit time due to the particles' 
departure from the site (loss flux) at site i, lf(t) is the 
gain flux, and nAi is the loss due to conversion. Particles' 
conservation for transitions between the two neighboring 
sites corresponds to 



Wi-i t ili-i(t) + w i+ i^I l+1 {t), 



(9) 



where Wij is a probability to jump to site j when leaving 

-M = = 1/2. 



i. For unbiased walks one has Wi 
Thus: 



^ = ^r-iW + i- it) K,Ai(t). (io) 

We now combine this continuity equation with the equa- 
tion for /" (t) following from the assumption about the 
distribution of sojourn times. The loss flux at time t is 
connected to the gain flux at the site in the past: the 
particles which leave the site i at time t either were at 
i from the very beginning (and survived without being 
converted into B), or arrived at i at some later time tf < t 
(and survived). A probability density that a particle 
making a step at time t arrived at its present position at 
time t' is given by the waiting time distribution ip(t — t'), 
the survival probability being p(t) — e~ Kt . Thus: 

lr(t)=Tf>(t)e- Kt M0)+ f '^(t-f)e- K ^lf(f)df. 

.In 

(11) 



Applying Eq.© we get: 
ir(t) = Mt)M0) 



dt' 



+ KAi{t') + ir {t') 



(12) 
dt' 



where ip s (t) = ip(t)e~ Kt is the non-proper waiting time 
density for the actually made new step provided the par- 
ticle survived. This approach can also be generalized to 
bimolecular reactions [ljj. Changing to the Laplace do- 
main and noting that ip s {u) = ip(u + k) we get 



I { (u) = <f> K (u)Ai(u) 



(13) 
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with given by 



(u + k) tp(u - 



1 — tp{u + k) 
Returning to the time-domain we thus get 



IT it) 



<$> K (t-t')Mt')dt', 



(14) 



(15) 



Note that $ K (t) given by the inverse Laplace transform 
of <l K (ii) corresponds to & K (t) = <& (t)e~ Kt where <I>o(i) 
obtained by taking k = is the usual memory kernel of 
the generalized master equation for CTRW. 

Combining Ea. lfT5|) with Eqs.© and © we get: 



dAj(t) 
dt 



*«(t - 



Ai-i(t') , A i+1 (f) 



- Ai(t') dt' - nAi(t). 
Transition to a continuum in space (x — ai) gives 



(16) 



dA(x,t) a 



dt 

,2 ft 



a 2 /"* 

— / $ K (i-<')AA(a;,i')*'-KA(x,t) 
2 Jo 



= tt / $o(* - Oe"' 5 ^''' A4(i, t')<it' - t), (17) 
2 Jo 

a rather unexpected form, where the reaction rate explic- 
itly affects the transport term. 

For the exponential waiting time distribution ip{t) = 
e~ xt corresponding to ip( u ) = A/(u + A) the kernel reads 
$o(£) = A5(t), and the existence of an additional expo- 
nential multiplier does not play any role: The reaction 
diffusion equation is perfectly exact. 

In the case of slowly decaying <I>o(i) the exponen- 
tial cutoff introduced by the reaction is crucial. For 
power-law waiting time distributions and for k = 
the integral operator J Q $o(i — t')f(t')dt' is the opera- 
tor of the fractional derivative: For such distributions 
tp(u) ~ 1 — (ru) a r(l — a) (where r is the appropri- 
ate time scale) and (for u — > 0) we have &o(u) ~ 
(1/t"T(1 — a))w 1-Q! which is proportional to the oper- 
ator of the Riemann-Liouville derivative of the order a: 



4 Jo $o(* - t')f(t')dt' = K a uD\- a f for sufficiently reg- 
ular functions /. The generalized diffusion coefficient 
reads K a = a 2 [2r a r(l — a)] -1 . For k > however the 
reaction affects the diffusion part of the equation: the 
Laplace transform of the integral kernel $ K (i) reads 



1 



T«r(l - a) 



(u + k) 



(18) 



and is no more a fractional derivative. The integral op- 
erator f t (l - a,K)f = T a T{l - a) Jo$ K (t - t')j{t')dt' 
corresponds in time domain to 



T t (l - a, k)S 



dt 



t g-K(t-t') 



t e - K {t-t') 



(19) 



turning to be a fractional derivative only for k = 0. The 
equation for the ^-concentration thus finally reads: 



dA(x,t) 
dt 



K a T t (l-a,K)AA(x,t) - nA(x,t). (20) 



Although our reaction does not depend on the particles' 
motion, the parameters of the reaction explicitly enter 
the transport operator T of the equation. 

Analogously we will now derive an equation for the B- 
particles. As for A one has a balance condition for the 
mean particle number Bi of B-particles on site i 



dBj(t) 
dt 



Jt(t)-Jr(t) + KMt), 



(21) 



where J 2 + denotes the gain flux and jf the loss flux of 
particles B at site i. The continuity equation reads: 



dBi{t) 



1 



A B-particle that leaves the site i at time t either has 
come there as a B-particle at some prior time or was 
converted from an ^4-particle that either was at site i 
from the very beginning or arrived there later, at t' > 0. 
Thus: 



Jr(t) = / ii(t-t')jf{t')dt' + iP(t)[l-e- Kt ]A i (0)+ / iP(t-t') 
Jo Jo 

lIBi( fL- KAi {t') + jr{t') 



1 - e 



-K(i-t') 



I+(t')dt' 



= / iit-t') 



+ / W-t') 

Jo 



dt 

I _ e -K(t-0 



,,Ai P- + KM t') + ir(t>) 



dt' 



dt' + iP(t) [l-e- Kt ] A t (0) + 
dt', 



(23) 
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where the local balance equation (JSJ was used. Applying Laplace transform and solving for J, (it) we get: 

J7(u) = ^— — \i)(u) uB t (u) - nAi(u) + -ijj(u) - tjj(u + k) uAi(u) + nAi(u) + 17 (u) X. 

1 — ip(u) L J L J L J > 



(24) 



Here the initial condition Bi(Q) 
Using Eq. l|13[) for 17 we get: 



was explicitly used. 



J7( u ) = $„(«)£<(«) + [$o(u) - (25) 
After inverse Laplace transformation we get: 



J7(t) = / $ (t-t')-B i (t / )rft / + 
Jo 

ft 



+ / *o(* - 







1 _ e -*(t-t') 



Ai(t')dt'. (26) 



We now substitute this into the continuity equation (|22H , 
perform the transition to a continuum and get 



dB(x,t) a 
di ~ 



2 rt 



$o (t - l/)AB(x, t')dt' + kA(x, t) + 



,2 r t 



1 - e 



-re(t-t') 



Ai(a:,t')A'. (27) 



For the exponential waiting time density for which 
$o(i) = \6(t) the third term in l|27|l vanishes and a nor- 
mal reaction-diffusion equation arises. For a power law 
waiting time distribution Eq. I|27(l can be written in terms 
of fractional derivatives and the operator T t (l — oi,k), 
Eq.lHU: 

dB(x, t) K ^ oD i-a AB ^ f j + KA(a , ; q + 



at 



£,i-«_f t (l-a,K) AA(i,t). (28) 



Note that the equation for the product contains the term 
depending on the concentration of the component A at 
all previous times. This term has to do with the fact 
that the products are introduced into the system later 
on in course of the reaction, and their motion therefore 
is described not by the normal CTRW (and the corre- 
sponding fractional diffusion equation) but by the aged 
one 0. Note also that the sum of Eqs.|J20J and 
always yelds the "normal" subdiffusion equation for the 
overall concentration C(x,t), Eq.©. Moreover the solu- 
tions of Eqs. and (gSJ satisfy Eq. ©. 

In summary, we derived here the equations describing 
the time evolution for the local A and ^-concentrations 
in a simple monomolecular conversion reaction A — > B 
taking place at a constant rate k and under subdiffusion 
conditions. These equations do not have a usual form 
of reaction-diffusion equations with the transport term 
independent on the reaction one. This fact is due to 



nonmarkovian property of the subdiffusion process, and 
will persist for more complex reaction schemes as well. 
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For example, for an A + B — > C reaction which 
takes place locally at a rate k the survival probabil- 
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ity p of an A particle in the time interval [t',t] is 
given by a solution of the equation dp/dt — —kB{€). 

It reads: p{t) = exp (- f* KBi(t')dt'^j . The loss 

flux of A-particles can then be written as follows: 

I- it) = V(t)exp(-/o«Bi(Od*')^(0) + />(* - 



t')exp(-f*KBi(t")dt")lf(t')dt'. Note that this ex- 
pression is a functional of the B-concentrations at all 
previous times. The overall structure of equations in this 
case is much more involved than in the monomolecular 
case considered in the present Letter. 



